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Simulations show that sliding bilayers of colloidal particles can exhibit a new phase, the "melt- 
freeze" phase, where the layers stochastically alternate between solidlike and liquidlike states. We 
introduce a mean field phenomenological model with two order parameters to understand the in- 
terplay of two adjacent layers while the system is in this remarkable phase. Predictions from our 
numerical simulations of a system in the melt-freeze phase include the tendency of two adjacent lay- 
ers to be in opposite states (solid and liquid) and the difference between the fluctuation of the order 
parameter in one layer while the other layer is in the same phase compared to the fluctuation while 
the other layer is in the opposite phase. We expect this behavior to be seen in future simulations 
and experiments. 

PACS numbers: 



I. INTRODUCTION 

General theories of nonequilibrium phenomena are still 
elusive and current research in nonequilibrium physics is 
directed on a phenomenon-by-phcnomenon basis. One 
important problem is the shear-induced melting (and/or 
freezing) of a solid. The shear flow of a solid is widely 
studied experimentally |lj due to its importance in ma- 
terial science. Being a standard example of the rich field 
of nonequilibrium phase transitions, shear-induced melt- 
ing is also widely studied theoretically on systems as di- 
verse as vortex lattices to suspended colloidal particles 

Mia. 

For clarity, we will focus on a specific system of col- 
loidal particles on two dimensional sheets 0. An iso- 
lated sheet at low enough temperature will have colloidal 
particles in an ordered phase. When two such sheets are 
driven across each other, one sheet subjects stress on the 
other (and vice versa) and at high enough stresses, could 
induce melting 4]. At large enough driving forces, the 
sheets are occasionally found to re-order to a re-entrant 
solid phase displaying, the phenomenon of shear-induced 
freezing [jj. 

Recently, numerical simulations of sliding bilayers by 
Das et al. discovered a new "melt-freeze" phase be- 

yond the standard shear-induced phases. Das et al. per- 
formed Brownian dynamics simulations of particles on 
two adjacent monolayers in the manner described in the 
previous paragraph. Earlier simulations usually contain 
one monolayer driven across a fixed substrate [4| which is 
applicable when one layer is much stiffer than the other. 
However the simulations of Das et al. addressed the sit- 
uation in which two comparably soft layers are driven 
across each other and thus the behavior of the particles on 
both layers becomes important. The relevant parameters 
in their simulations included the strength of the driving 
force, the magnitude of the coupling between particles 
on differing sheets, and the noise amplitude. Changing 
these parameters displays rich and unexpected behavior. 

Das et al. found that at fixed interlayer coupling and 
noise strength, the sytem undergoes different phases as 



the driving force is varied. For small driving forces, the 
two monolayers remain in the ordered phase and simply 
creep past one other. The same qualitative behavior is 
seen for very large driving forces. However, at interme- 
diate driving forces, the two monolayers undergo a melt- 
freeze phase in which both layers stochastically disorder 
(melt) and order (freeze). As the interlayer coupling in- 
creases, the melt-freeze phase persists over a larger range 
of driving forces. For a fixed driving force (in which the 
system is in the melt-frcczc phase), increasing the cou- 
pling also increases the amount of time a layer spends in 
the disordered state. 

In addition to the simulation of particles on layers, 
Das et al. introduced a simple mean field phenomenologi- 
cal model . This phenomenological model contains one 
order parameter and a strain variable. Like the spatially 
dependent system, the mean field model has three simu- 
lation parameters (noise strength, coupling strength, and 
driving force). Remarkably, this simple model contains 
much of the rich behavior of the more detailed parti- 
cle simulations. However, the interplay between the two 
layers cannot be studied by this model with only one or- 
der parameter. The simple model also clearly connects 
the melt-freeze phase with stochastic resonance phenom- 
ena H0. 

In this paper^we closely follow but extend the phe- 
nomenology in [4| to include two order parameters p\ 
and pi in the same spirit that Das et al. 's particle simu- 
lations 0, @ generalized the conventional models of one 
sheet sliding over a fixed substrate. This natural exten- 
sion qualitatively recovers much of the behavior seen in 
the one order parameter phenomenological model (there- 
fore also in the particle model). In addition, we observe 
two new features involving how one order parameter af- 
fects the free energy landscape of the other. The first 
feature is the tendency of the two order parameters to be 
in different phases when the coupling is finite, i.e. when 
pi is in the ordered phase, pi is likely to be in the disor- 
dered phase. The second feature involves the fluctuation 
of the order parameters. When p\ and p2 are out of 
phase, their fluctuations are smaller than when they are 
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in the same phase. These features are compared to the 
simulation of particles in 0, [fj . 

In the next section, we describe our two parameter 
Landau model and the Langevin-Ginzburg dynamics. 
Section III contains the results from our simulations. As 
in 0,0, the results include driving force dependence and 
coupling strength dependence. 

II. MODEL FREE ENERGY AND DYNAMICS 

The order parameter in our mean field model repre- 
sents the structure factor for the colloidal particles in a 
given sheet The model free energy must represent a 
system with two phases (crystal and liquid) and contain 
a first order transition between the phases. A simple ef- 
fective Landau free energy with these considerations is 
the simple double well potential used in Ref . 0, ||| , 

V(p) = f P 2 - f p 3 + \ P \ (1) 

where a, /?, and 7 are all positive. The values are chosen 
so that there are two phases (two minima) and the stable 
phase is the crystalline one when there is no stress. 

Modelling the dynamically driven phase transitions re- 
quires a strain variable 6, as in 0, |(|. Since stress is 
induced by a periodic crystal, the strain variable must 
be periodic (we use the range [0, 1)). A liquid does not 
induce stress on a driven solid sheet. Using these proper- 
ties of strain, we model the part of the free energy with 
the strain variable to be 

W( P i,p2,ff) = lf%f%[l-cos(2*6)], (2) 

where c is a coupling strength. The total model free 
energy includes both order parameters in a free energy 
landscape of Eq. Q and the strain variable and coupling 
ofEq.©, 

F(p 1 ,p2,6)=V(p 1 ) + V(p 2 )+W( Pl ,p 2 ,6). (3) 

Figure n shows the free energy for an order parameter 
while the other is in the ordered state. Since the other 
sheet is ordered, it can induce stress. As stress increases, 
shown by a larger 9, the ordered phase becomes less sta- 
ble and the system is driven towards the liquid phase. 
The changing landscape shown in figure^shows how the 
simple free energy of Eq. @ can represent shear- induced 
melting. 

Eq. displays the equilibrium conditions in which the 
order parameters exist. For these driven systems, the dy- 
namics are also important. In this paper, we use over- 
damped Langcvin dynamics of the form 

1 dF 

Pl = ~ p dp^ + ^ 
1 dF 

r p dpi 

1 dF 

= -^^E+ d + Ve- (4) 
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FIG. 1: The free energy for one order parameter p\ while the 
other order parameter pi is in the ordered state. The different 
plots represent different values of 9 from 0.0 (no stress) to 
0.5 (maximal stress). The coupling c is set to 0.25 here and 
a = l,/3 = 6.25, and 7 = 8. 



The damping coefficients T's are set to unity, d repre- 
sents a constant driving force for the strain, and the 77's 
represent noise. We use Gaussian noise with zero mean 
< 77 >= and variance D, < r){t)r)(t') >= 2DS(t - t'). 
The amount of noise D is related to the temperature 
through the fluctuation-dissipation theorem. The strain 
variable and the order parameters have noise of equal 
amplitudes in our simulations. 

Eq.© and Eqs.QJ describe the free energy and dy- 
namics of our phenomenological model. In the next sec- 
tion, we briefly describe the numerical implementation of 
these equations and the results of our simulations. 



III. RESULTS 

We numerically integrate the Langevin-Ginzburg equa- 
tions, Eqs.QJ, using the fourth order Runge-Kutta algo- 
rithm. For all simulations, the stepsize is 0.001 and the 
integrator is run for ~ 10 6 transient steps which are dis- 
carded. We also set a = X,j3 = 6.25, and 7 = 8 for all 
simulations. As in 0,0], the adjustable parameters in 
the model are the strength of the constant driving force 
d, the coupling strength c, and the noise amplitude D. 
Below, we show how the system behaves by varying these 
parameters independently. To study the variation of the 
driving force and the coupling, we keep the noise am- 
plitude fixed at D = 0.004. We observe that this noise 
strength is insufficient to drive the system to the disor- 
dered phase from the ordered phase for very long times 
when stress is absent. 
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A. Driving Force Dependence 



As described in the introduction and in |J, |6j , the sys- 
tem stays essentially crystalline when the driving force is 
below a threshold value or above another threshold. Only 
at intermediate driving forces is the melt-freeze phase 
seen. In figure El we show pi and p2 versus time at con- 
stant coupling, c = 0.25. The unit of time in figure Eland 
figure 03 represents 2 x 10 4 timesteps. Figure |3 shows the 
distribution of p\ (the distribution of pi is the same after 
averaging over long enough time due to the symmetry 
between the two in the free energy). 

To understand this behavior qualitatively, we discuss 
the role of the strain variable in determining the phase of 
the system. The driving force tilts the cosine potential 
and drives the strain variable in a preferred direction. 
At very small driving forces, the strain variable is rarely 
able to leave the minimum in which it started. As seen 
from figure ^ when 9 is near this minimum the crystal 
well remains stable and p does not cross the barrier. For 
very large driving forces the cosine potential is not seen 
and the distribution of 9 is essentially flat and does not 
spend enough time near its peak to allow p to cross over 
to the liquid well. In between these two limits exists an 
optimal value for d. Figure 0] shows the distribution of 9 
(since 9 is periodic from zero to one, it is only necessary 
to display the distribution of the decimal) for different 
driving forces. As d increases, the peak in the distribu- 
tion moves to a greater value and the whole distribution 
flattens. 

From these figures, it is clear that there is an optimal 
driving force in which the melt-freeze phase exists. This 
behavior qualitatively agrees with both the particle sim- 
ulations and_the one order parameter phenomenological 
model in 



B. Coupling Dependence 

Keeping the driving force at a fixed value of d = 0.1 
(which is near the optimal value), we observe the behav- 
ior of the system as the coupling constant varies. Fig- 
ure[H]sh.ows the order parameters versus time and figure|S| 
shows the distribution of p\ for different couplings. 

Since the driving force is taken to be near the optimal 
value, all three cases from figures [S] and are in the melt- 
freeze phase. From these figures, it is clear that as the 
coupling strength increases, the time spent in the liquid 
phase increases as well, which is a feature seen in the 
models of Das et al. Furthermore, systems with larger 
coupling strengths have order parameters that tend to 
be in opposite phases. 

To explore this further, we define the boundary be- 
tween the two phases to be the location of the local max- 
imum without stress, p = (l/2-y)[f3 ~ (f3 2 - -kry) 1 / 2 ], and 
calculate how often p\ and pi are in the same well. In 
the weak coupling case (c = 0.2), the two order parame- 
ters are in the same well approximately 80 percent of the 





FIG. 2: The order parameter vs. time for three different 
driving forces d = 0.025, 0.1, and 10.0 respectively from top to 
bottom. The coupling c = 0.25 and noise strength D — 0.004 
are kept fixed. 



time. For the medium coupling case (c = 0.25), p\ and 
P2 are in the same phase approximately 35 percent of the 
time and in the large coupling case (c = 0.3), they only 
reside in the same well about 10 percent of the time. In 
the infinite coupling limit, the two order parameters will 
"repel" each other rather strongly. Due to the symmetry 
in the free energy between p\ and p%, in this limit, each 
order parameter has the same probability of being in the 
ordered phase as the disordered. 

The observation described in the above paragraph pro- 
vides a new prediction for spatially dependent Simula- 
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FIG. 3: The distribution of p\ taken over 10 1 steps for d — 
0.025, 0.1, and 10.0. The coupling c = 0.25 and noise strength 
D = 0.004 are kept fixed. 
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FIG. 4: The distribution of 8 taken over 10 11 steps for d — 
0.025,0.05,0.1, and 10.0. The coupling c = 0.25 and noise 
strength D — 0.004 are kept fixed. 8 is periodic in the range 
[0,1). 





5000 10000 
lime 



15000 20000 




tions and experiments. Obviously, the behavior seen here 
cannot be seen in the one order parameter mean field 
model. 

By observing that the coupling term in the free energy 
merely changes the coefficient of the quadratic term for 
a given order parameter, one can intuitively understand 
this behavior. In the large coupling regime, the landscape 
in which p\ exists is strongly determined by pi and 8. If 
both start in the crystal well, as soon as 8 is finite, this 
minimum quickly becomes unstable and one of the p's 
must roll down to the disordered minimum. If p2 falls 
down the liquid well, the amplitude of the effective cou- 
pling drops near zero and p\ is again happily in its stable 
well. In addition, from p2 > s perspective, it is also in its 
own stable well because cp\/2 is finite and large. Since 
the height of the cosine potential is proportional to both 



FIG. 5: The order parameter vs. time for three different 
coupling strengths c = 0.2, 0.25, and 0.3 respectively from 
top to bottom. The driving force d = 0.1 and noise strength 
D = 0.004 are kept fixed. 



order parameters squared, when one order parameter is 
in the liquid phase, 8 is essentially driven without resis- 
tance. Driving the strain variable in this way allows the 
crystal well to become stable for p^. With aid of noise, 
P2 can eventually move back to the crystal phase. But 
once this occurs, both order parameters are in a highly 
unstable position and one will quickly leave for the liquid 
phase and the situation repeats itself. Which one crosses 
to the liquid phase first is a stochastic process and cannot 
be determined, but the repulsion effect still remains. 
The simulations of layers of particles in 0, [|| show 
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FIG. 6: The distribution of pi taken over 10 11 steps for c = 
0.2, 0.25, and 0.3. The driving force d = 0.1 and noise strength 
D = 0.004 are kept fixed. 



hints of this repulsive behavior though this repulsion is 
not as explicitly seen as in our mean field model. Our 
phenomenological model provides intuition into this be- 
havior and predicts that it should be more easily seen as 
the coupling increases. More work is necessary to verify 
whether this is seen in a spatially dependent model. 

With the form of our model free energy and the pa- 
rameters in our simulations, the repulsive behavior is 
so strong that the system has yet to be seen in a state 
where both order parameters are in the disordered phase. 
Though this has not been seen in our simulations, the 
stochastic nature allows this state to be possible and, 
presumably, parameters can be tuned so that both or- 
der parameters can be in the liquid well in a quasi-stable 
state. 



C. Fluctuations 

A close inspection of figures [21 and [S] reveals interesting 
behavior regarding the fluctuations of the order parame- 
ters; the fluctuation of one order parameter depends on 
the state of the other order parameter. Obviously, the 
noise amplitude and shape of the well determines the 
amount an order parameter fluctuates. The magnitude 
of the other order parameter (and the strain) influences 
the quadratic term of the free energy, thereby changing 
the shape of the well. 

Figures |21 and show that the fluctuations of an order 
parameter in the crystalline well is much larger when 
the other order parameter resides in the same well. To 
quantify this, we compare the standard deviation of p\ 
while both p\ and p2 are in the crystalline well (p > 
(I/27) \J3- (/3 2 - 4aj) 1/2 ] ) to the standard deviation while 
pi is ordered and pi is disordered. For c = 0.25, d = 0.1, 
and D — 0.004, the standard deviation is approximately 
0.039 for the former and 0.028 for the latter. 



This behavior makes physical sense because when both 
layers are solid, stress is induced between the two layers 
and the particles are driven further away from their equi- 
librium positions. Whereas, if the particles in one layer 
are disordered, the layer will not induce any stress on the 
neighboring layer. 



IV. DISCUSSION 

We note two minor disagreements between our mean 
field model and the more detailed particle simulations 
of Das et al. (3, El- The first is that Das et al. ob- 
served a smaller extent of order with smaller interlayer 
coupling whereas we observe the opposite, i.e. for our 
model, the mean value for p in the crystalline state is 
larger for smaller coupling c. The second disagreement 
is that for larger coupling we see larger fluctuations of 
the order parameter while Das et al. observed a decrease 
in the amount of the fluctuations. Considering the sim- 
plicity of the phenomenological model, disparities such 
as these are expected and the many agreements are quite 
surprising. 

By varying the noise in the one order parameter phe- 
nomenological model, Das et al. connects this melt-freeze 
phase to stochastic resonance |4]. Qualitatively, the op- 
timal driving force exists when the driving time constant 
matches (within perhaps a factor of two) the thermal 
barrier crossing time constant. Conversely, with a fixed 
driving force, one expects a small window of noise am- 
plitudes where the thermal escape rate "resonates" with 
the driving rate. The existence of an optimal noise am- 
plitude 0| is a signature of stochastic resonance Q . Due 
to the similarities of that model and ours, we expect that 
our model displays stochastic resonance as well, though 
we have not performed extensive simulations with noise 
dependence. 

In this paper, we introduce a simple phenomenological 
model for a system of two sheared monolayers to study a 
novel phase, the melt-freeze phase, first seen in the simu- 
lations of Das et al. 0, 01 • Using numerical simulations, 
we observe that the two order parameter model recovers 
many of the same qualitative features as the simulations 
from Das et al. , especially the dependence on driving 
force and coupling strength. Our simple model displays 
new behavior involving the interplay of the order param- 
eters of two adjacent layers. One behavior is that with 
strong interlayer couplings, two adjacent layers will be in 
opposite states of matter. Another involves the fluctu- 
ation of the order parameters in the melt-freeze phase; 
fluctuations are larger when both order parameters co- 
exist in the same well than when they are in opposite 
wells. More simulations on spatially dependent systems 
are necessary to verify these observations. 

Das et al. suggested experiments using bulk colloidal 
crystals, ordered copolymer monolayers, or colloidal 
monolayers to see this unusual melt-freeze phase |4| . Sys- 
tems of charged colloidal particles are good candidates 
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to observe this phenomenon because they possess many 
properties that are tunable. The review article by Pal- 
berg thoroughly discusses the physical properties 
of charged colloids, how to change the properties, and 
different measurement techniques. For instance, the in- 
terparticle interaction can be adjusted by changing the 
electrolyte concentration, thereby altering the Coulomb 
screening. As seen in the simulations, the melt-freeze 
phenomenon exists only in certain parameter windows, 
therefore tunable parameters are necessary for experi- 
mental observation. The melt-freeze phase and the be- 



havior seen in our mean field model could be seen in these 
experimental systems. 
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